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Abstract 

Self-organization and pattern formation in network-organized systems emerges from the eolleetive 
aetivation and interaetion of many intereonneeted units. A striking feature of these non-equilibrium 
struetures is that they are often loealized and robust: only a small subset of the nodes, or eell assembly, 
is aetivated. Understanding the role of eell assemblies as basie funetional units in neural networks and 
soeio-teehnieal systems emerges as a fundamental ehallenge in network theory. A key open question is 
how these elementary building bloeks emerge, and how they operate, linking strueture and funetion in 
eomplex networks. Here we show that a network analogue of the Swift-Hohenberg eontinuum model— 
a minimal-ingredients model of nodal aetivation and interaetion within a eomplex network—is able to 
produee a eomplex suite of loealized patterns. Henee, the spontaneous formation of robust operational 
eell assemblies in eomplex networks ean be explained as the result of self-organization, even in the 
absenee of synaptie reinforeements. Our results show that these self-organized, loeal struetures ean 
provide robust funetional units to understand natural and soeio-teehnieal network-organized proeesses. 
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INTRODUCTION 


Pattern formation in reaction-diffusion systems [HI |3 has emerged as a mathematical 
paradigm to understand the connection between pattern and process in natural and sociotech- 
nical systems (Si. The basic mechanisms of pattern formation by local self-activation and lat¬ 
eral inhibition, or short-range positive feedback and long-range negative feedback [|4l El are 
ubiquitous in ecological and biological spatial systems, from morphogenesis and developmen¬ 
tal biology [HI [6l to adaptive strategies in living organisms [[71 [H and spatial heterogeneity in 
predator-prey systems 0. Heterogeneity and patchiness in vegetation dynamics, associated 
with Turing patterns in vegetation dyanmics have been proposed as a connection between pat¬ 
tern and process in ecosystems [HOllIIl, suggesting a link between spatial vegetation patterns 
and vulnerability to catastrophic shifts in water-stressed ecosystems iTTTl - fTdl . 

The theory of non-equilibrium self-organization and Turing patterns has been recently ex¬ 
tended to network-organized natural and socio-technical systems [HEHIHl, including complex 
topological structures such as multiplex [HU |20l, directed [|2T1 and cartesian product net¬ 
works [[22l . Self-organization is rapidly emerging as a central paradigm to understand neu¬ 
ral computation Il23l - E5[| . The dynamics of neuron activation, and the emergence of collective 
processing and activation in the brain, are often conceptualized as dynamical processes in net¬ 
work theory [[26l - l28l . Self-organized activation has been shown to emerge spontaneously from 
the heterogenous interaction among neurons [[25ll . and is often described as pattern formation 
in two-population networks [[29l - [32ll . Localization of neural activation patterns is a conceptu¬ 
ally challenging feature in neuroscience. Cell assemblies, or small subsets of neurons that fire 
synchronously, are the functional unit of the cerebral cortex in the Hebbian theory of mental 
representation and learning [[33I4371 . Associative learning forms the basis of our current un¬ 
derstanding of the structure and function of neural systems [[38l - l40ll . It is also the modeling 
paradigm for information-processing artificial neural networks ll4TI - l43ll . The emergence of cell 
assemblies in complex neural networks is a fascinating example of pattern formation arising 
from the collective dynamics of interconnected units Il25ll4^ . Understanding the mechanisms 
leading to pattern localization remains a long-standing problem in neuroscience [[341.144l - l48l . 

Here we show that simple mechanisms of nodal interaction in heterogeneous networks allow 
for the emergence of robust local activation patterns through self-organization. The simplicity 
and robustness of the proposed single-species pattern-forming mechanisms suggest that analo- 
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gous dynamics may explain localized patterns of activity emerging in many network-organized 
natural and socio-technical systems. We demonstrate that robust local, quantized activation 
structures emerge in the dynamics of network-organized systems, even for relatively simple 
dynamics. We propose a minimal-ingredients, phenomenological model of nodal excitation 
and interaction within a network with heterogeneous connectivity. Our goal is to demonstrate 
that a simple combination of local excitation of individual units, combined with generic exci¬ 
tatory/inhibitory interactions between connected units, leads to self-organization, and can ex¬ 
plain the spontaneous formation of cell assemblies without the need for synaptic plasticity or 
reinforcement. Our model can be understood as a network analogue of the Swift-Hohenberg 
continuum model [|49ll49l - [5TI . and is able to produce a complex suite of localized patterns. The 
requirements are minimal and general: simple local dynamics based on canonical activation po¬ 
tentials, and interactions between nodes that induce short-range anti correlation and long-range 
correlation in activation. Because of their robustness and localization, self-organized structures 
may provide an encoding mechanism for information processing and computation in neural 
networks. 

MODEL OF NETWORK DYNAMICS AND STABILITY ANALYSIS 

We restrict our analysis to the simplified case of symmetric networks, but our main results 
can be generalized to other network topologies, including directed ll^ and multiplex [|T^ l20l 
networks. A node’s state of activation, measured through a potential-like variable u, is driven by 
local excitation dynamics and by the interaction with other nodes in the network via exchanges 
through the links connecting them. In dimensionless quantities, the proposed excitation- 
inhibition model for the evolution of potential, Uj, in each node i = 1,..., A^, is given by 
the model 

^ = /(Mi,/i) + Ji, (1) 

where /(«*, /i) is a dynamic forcing term, representing a double well potential, and /r is a bi¬ 
furcation parameter that will be used to establish the conditions for stability and localization of 
the response patterns (Fig. 1 A). The currents, Jj, represent the excitatory/inhibitory interactions 
among nodes in the network. The structure of these nodal interactions is one of the key pat¬ 
tern forming mechanisms in the present model. We consider short-range anti-correlation, and 
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higher-order, longer-range dissipative interaetions. This two-level interaetion strueture, whieh 
induees anti eorrelation in the short range (nearest-neighbors, or first-order eonneetivity), and 
long-range eorrelation (seeond-nearest neighbors, or seeond-order eonneetivity) is represented 
in Fig. lA. Mathematieally, we express the integration of synaptie eontributions as 

N N 

= ( 2 ) 

i=i i=i 

The strueture of the above nodal interaetions turns the dynamies ([T])-@ into a network anologue 
of the Swift-Hohenberg eontinuum model, 

f)ll 

— = /K/i)-(l + VY«, (3) 

whieh is a paradigm for pattern-forming systems [l49l - [5^ . The simplest form for the in¬ 
teraetion matriees representing these eorrelation/anti-eorrelation effeets (while ensuring that 
the interaetion fluxes eonserve mass or eharge) is based on network representation of Lapla- 
eian and bi-Laplaeian operators, and respeetively. The network Laplaeian, 
is a real, symmetrie and negative semi-definite NxN matrix, whose elements are given 
by = Aij — kiSij [[HI, where Aij is the adjaeeney matrix of the network, ki = 
the degree (eonneetivity) of node i and 5ij is the Kroneeker delta. A diffusive, Fiekian-type flux 
of the aetivation potential u to node i is expressed as Yl!j=i ^ ~ 

Figure [T]4-top). Plain waves and wavenumbers on a network topology are represented by the 
eigenveetors <p°‘ = .., 0^^^) and the eigenvalues Aq, and of the Laplaeian matrix, whieh 

are determined by the equation with a = 1,.., [[2l. All eigenvalues 

are real and non-positive and the eigenveetors are orthonormalized as 
where a, [3 = 1,..., N. The elements of the bi-Laplaeian matrix of a network ean be expressed 
as Llf = YlZi where the {A^)ij = Hz AiAj matrix 

has information about seeond order nodal eonneetivity and takes nonzero values if node i is two 
jumps away from node j. The operation models negative diffusion (inhibition) 

from the first neighbors of node i and at the same time diffusion from its two-jump neighbor¬ 
hood (see Figure [^-bottom). The bi-Laplaeian, has the same eigenveetors as (i.e. 
0") and its eigenvalues are the square of those of L^'^\ A^. 

To understand the properties and pattern-forming meehanisms in our model, we first inves- 
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tigate the stability of flat states of the dynamieal system 


dui 

dt 


N 


N 


f{ui, /i) - 2 ^ Lf^Uj -^L 


( 4 ),, 


i=i 


i=i 


i = 1, 


(4) 


Flat, stationary solutions u of Eq. (4) satisfy f{ui, jj) = 0, where the nodal state of aetivation 
is equal for all nodes in the network, Ui = u,\/i = 1, For f{ui,^) = —(1 + ^)ui + 

1.5u^ — uf, there are three uniform solution branehes given by mq = 0 and u± = [1.5 ± 
■y/ 1.52 — 4(/i + l)]/2. It is well known in one and two dimensional continuum spaces that 
these uniform states can become unstable and a wealth of self-organized patterns can arise [|49l - 
[5^ . In a linear stability analysis, the stability of flat stationary solutions to small perturbations 
is determined by the eigenvalues of the Laplacian and bi-Laplacian matrices. Introducing small 
perturbations, 6ui, to the uniform state u, Ui = u + 6ui, the linearized version of Eq. (4) 
takes the form ddui/dt = fu{u)6ui - where /„ = duf = 

. After expanding the perturbation 5ui over the set of the Eaplacian eigenvectors, 5ui = 
Z]q=i where is the expansion coefficient, the linearized equation is transformed 

into a set of N independent linear equations for the different normal modes: 


Aa = fu{u) - 2A„ - A^, a = 1,..., A^, 


(5) 


where Aq, are the eigenvalues of the Eaplacian matrix. The a-mode is unstable when Re Aq, is 
positive. Instability occurs when one of the modes (the critical mode) begins to grow. At the 
instability threshold. Re Aq, = 0 for some Oc and Re A < 0 for all other modes. In Figure IB 
and C we summarize the linear stability analysis of the flat states of our model on a scale-free 
network constructed using the Barabasi-Albert model (BA) of network growth and preferential 
attachment [f54l . We find that, indeed, there is a large parameter range for which the resting 
potential is stable. As we demonstrate below, in the stable regime, input stimuli may trigger 
localized patterns of activation. 


LOCALIZED PATTERNS 

Localized activation patterns are possible due to the particular structure of the model, with 
short- and long-range nodal interactions. Mathematically, the localized states are homoclinic 
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orbits in the network spaee around the base resting state, u = uq = 0. The existenee of these 
homoelinie orbits ean be studied using the teehnology developed for the linear stability analysis. 
Sinee homoelinie orbits leave the flat state as we approaeh a small neighborhood (eluster) of 
the network, the fixed point must have both stable and unstable eigenvalues. We linearize 
Eq. (4) around u = 0 and write 5ui = Qa ^ 1, arriving to the relation /„(0) — 

2Aq — = 0. Sinee the Laplaeian eigenvalues A^ are real and non-positive values, we ean 

write them in the form Aq, = — If /i > 0 the topologieal eigenvalues of u = 0 form a 
eomplex quartet ka = ± ^ -f C>(/i). For /i = 0 they eollide pairwise on the imaginary 

axis, and for /i < 0 they split and remain on the imaginary axis ka = ^(±1 ± + 0{^). 

For /i = — 1 two of the topologieal eigenvalues eollide at the origin and for /x < — 1 they move 
onto the real axis. These results are summarized in Figure]^. The topologieal eigenvalues in 
the neighborhood of /i = 0 are eharaeteristies of the reversible 1 : 1 resonanee bifureation. 
Theory shows that under eertain eonditions the hyperbolie regime eontains a large variety of 
topologieally loealized states llSTl . 

To understand the onset of loealized patterns for different model parameters and input stim¬ 
uli, we eonstruet the bifureation diagram of the resting state, as a funetion of the total potential 
energy of the stimulus and bifureation parameter /x, in the vieinity of /x ~ 0 llSTl [52ll . A sin¬ 
gle bifureation braneh—eonstrueted using a pseudo-arelength eontinuation method f[55l —has a 
eharaeteristie “snaking” strueture of loealized states with varying aetivation energy ||u°||i 2 = 

.[^). As the system jumps from one steady state braneh to the next one, 
a new neighborhood in the network is being aetivated. Figurevisualizes the different steady 
loealized states of the six different branehes as they are spotted in the diagram of Figure [^. 
The response of the system is quantized: the transition from one pattern of aetivation to another 
one is diseontinuous as we vary the aetivation energy | |u°| |, or the parameter /x (Fig.|^). These 
jumps in aetivation energy eorrespond to the addition of neighbor nodes to the eluster (Fig.|^). 

The diseontinuous—quantized—nature of the network response leads to robustness in the 
loeal, final equilibrium patterns with respeet to the input signal amplitude. To gain insight into 
the robustness of the loealized patterns of aetivation, we performed a synthetie test in whieh 
we initially stimulate a speeifie neighborhood in the network, where we set xx* = ?x > 0 (i.e. a 
step-like funetion signal in network topology) and let the system evolve to equilibrium without 
deeay. We gradually inerease the amplitude u of the initial signal, and reeord the final energy 
values of the equilibrium, loealized states. For small amplitudes the perturbation relaxes baek 
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to the resting state, and no activation pattern is elicited. There is a threshold in the energy of 
the input stimulus beyond which robust quantized states are form. The states are robust in the 
sense that further increments in the input signal amplitude do not change the final equilibrium 
pattern (Fig.[^). 

The self-organized local structures are also robust with respect to random noise in the initial 
stimulus. We perform Monte Carlo simulations that probe the impact of the noise-to-signal ratio 
on the energy of the emerging quantized state. We have confirmed that the presence of small- 
amplitude noise has no effect on the equilibrium states of nodal activity. As can be expected, we 
do observe a departure from the energy of the base equilibrium state when the noise-to-signal 
ratio is sufficiently large, thereby masking the base stimulus altogether (Fig.[^). 


MEAN-FIELD APPROXIMATION OF THE GLOBAL ACTIVATION PATTERNS 


Our model predicts a range of parameter values where localized states disappear, and are 
replaced by global activation patterns. Mathematically, global patterns are possible when the 
non-active stationary solution is perturbed outside the parameter region of localized patterns 
(/i < 0). These—global—Turing patterns [[B HH can be understood and modeled using the 
Mean-Field Approximation (MFA), a method that segregates nodes according to their degree 
and has been successfully used to approximate a wide variety of dynamical processes in hetero¬ 
geneous networks, like epidemic spreading [|56l - [5^ . activator-inhibitor models [fT^ and voter 
models [[59ll . 

This theory allow us to reduce the problem to a single equation for the membrane potential 
for all the nodes in the system. Since in our model both the degree and two-jump degree play 
important role in the formation of patterns, we use a MFA where we assume that all the nodes 
with the same degree and two-jump degree behave in the same way. We start by writing Eq. (3) 
in the form 


dui 

dt 


f {li Qi kihi “f 


( 6 ) 


where the local fields felt by each node, hi = 9i = 

^ij{kjUj) are introduced. These local fields are then approximated as h* ~ k ~ 

k^Huu and gi ~ where ki = is the degree and /cf^ = is 

the number of secondary connections of node i (two-jump degree). The global mean fields 
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are defined by Hu = il/N) Y,k ^kHu'^ where = {l/{kNk)) Y.iek Y.j and Huu = 
(1/^) Ea:( 2 ) where Hif'^ = N,m)) Here, is the 

number of nodes with degree k, Ni^(2) is the number of nodes with number of two-jump 
neighbors and N = = Efc(2) ^k(^) the size of the network. In the above expressions, 

with Eiefc denote the sum over the nodes with degree k and with Eiefc(2) the sum over the 
nodes with two-jump nodal eonneetivity k"^. 

With this approximation, the individual model equation on eaeh node interacts only with the 
global mean fields Hu and Huu and its dynamics is described by: 

(111 

- = f{u) -2a{Hu -u)- (7) 

- [/3Huu-a‘^Hu-/3Hu +a'^u]. 

Since all nodes obey the same equation, we have dropped the index i and introduced the param¬ 
eters a{i) = ki and j3{i) = kf\ The activation potential depends now on the global fields Hu 
and Huu as well as on the parameter compination {a, (5), i.e. u = u{Hu, Huu, o:, l3). If the 
global mean fields Hu and Huu are given, the combination (a, [3) plays the role of a bifurca¬ 
tion parameter that controls the dynamics of each node in the system. The time independent 
version of above mean field equation can be written as a third degree algebraic equation that 
we solve N times for the N nodes in the system. For each node i, we get three solutions u\, 
/ = 1..3 that can be stable or unstable depending on the sign (negative or positive respectively) 
of the operator f'\ui +2a — a^. 

After tuning the bifurcation parameter /i to a negative value, we can compute the global 
Turing pattern from direct numerical simulations and determined the global mean fields Hu 
and Huu- Each node i in the network is characterized by its degree and second nodal connec¬ 
tivity, so that it possesses a certain parameter combination, (a, [3). Substituting these computed 
global mean fields as well as the values of a and [3 into equation (7), bifurcation diagrams of a 
single node can be obtained and projected onto the Turing pattern. In Figure (4A) we show for 
our “toy network model” that the stable brunches of the nodal bifurcation diagrams calculated 
using the MFA fit very well the computed Turing pattern. We further assess the dependence 
of the network topology on the global pattern formation and we find that when the degree dis¬ 
tribution is narrower compared to a scale-free network, the distribution of the (q;,/ 9) is more 
homogeneous and therefore the stationary Turing patterns look smoother (Figure 4B and C). 
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Therefore, global network Turing patterns are essentially explained by the bifurcation diagrams 
of individual nodes coupled to the global mean fields, with the coupling strength determined by 
their degree and two-jump connectivity. 


CONCLUSIONS 

Our results suggest a new mechanism for the formation of localized nodal assemblies in 
networks, arising from long-range—second neighbor—interactions. Rather than relying on re¬ 
inforcing mechanisms—synaptic plasticity, we show that localized, robust nodal assemblies 
are possible due to self-organization. The emergence of localized activation patterns derived 
from the simple and general functional structure of our proposed conceptual model: local dy¬ 
namics based on activation potentials, and interactions between nodes that induce short-range 
anticorrelation and long-range correlation in node-to-node exchanges. The proposed system is 
a network analogue of the Swift-Hohenberg continuum model, and is able to produce a complex 
suite of robust, localized patterns. Hence, the spontaneous formation of robust operational cell 
assemblies in complex networks can be explained as the result of self-organization, even in the 
absence of synaptic reinforcements. Hence, these self-organized, local structures can provide 
robust functional units to understand natural and technical network-organized processes. 
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Eellowship from the Spanish Ministry of Economy and Competitiveness (to ECE). 
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FIG. 1: (A) Pictorial illustration of our dynamical model of network interactions. Locally, the nodal acti¬ 
vation state is driven by the dynamic forcing term f{u, fi). In the inset we show the potential landscape— 
minus the integral of / with respect to u, which exhibits a single well (at u = 0) with an inflection point, 
a necessary condition for localized patterns to exist. Nodes interact in the network through diffusion-like 
exchanges via the links connecting them. The network Laplacian operator, represents short range 
diffusion of the species in the system (top). The network bi-Laplacian operator, induces short range 
anti-correlation with the nearest-neighbors, and long-range correlation with the second-nearest neighbors 
(bottom). (B-C) Linear stability analysis of the flat stationary solutions of our model. (B) The maxi¬ 
mum value of the growth rate A as a function of the bifurcation parameter fi for the two flat stationary 
states (brown) and u- (blue) on a Barabasi-Albert network model with mean degree {k) = 3 and 
size A^=2000. When the maximum value of A is negative, the state is stable with respect to small non 
uniform perturbation. (Inset) The growth rate A as a function of the Laplacian eigenvalue A (Eq. (5)) 
for three different values of the bifurcation parameter /r as they indicated in the main diagram for the 
flat stationary solution u_. (C) The flat stationary solutions uq and u± as a function of /x on the same 
network. Solid (dotted) lines represent stability (instability) with respect to small non-uniform perturba¬ 
tions. The labelled bifurcation points are /xq = 0, /xi = —0.44 and /x+ = —0.62 and /x_ = —1.82. The 
pink shaded region is where we observe localized self-organization patterns with respect to the trivial 
solution Mo- For values of ^ outside that region we get either global activation patterns (for /x < /xq) or 
any perturbation relaxes back to the flat stationary solution (for /x > fi*). 
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FIG. 2: Localized self-organized quantized patterns. (A) Stability of the trivial flat stationary state of 
our model with respect to the values of the bifurcation parameter, p,. For positive values of p the trivial 
stationary solution is stable with respect to uniform small random perturbations (solid line) while for 
negative values of p this state becomes unstable (dotted line). Also shown in the insets are the topological 
eigenvalues of the trivial state as we tune the bifurcation parameter. The behavior of eigenvalues in the 
neighborhood of /r = 0 indicates the possibility for localized patterns in the neighborhood of small 
positive values of p (pink shaded region). (B) A single branch of the bifurcation diagram in a Barabasi- 
Albert network model of size N=200 with mean degree equal to {k) = 3 and minimum degree equal 
to 1. Solid (dotted) lines represent stable (unstable) localized solutions. (C) Visualization of the localized 
patterns corresponding to the states indicated on the bifurcation diagram (B). Gray-colored nodes are 
non-active (rt = 0), red-colored nodes are active with rt > 0 and blue-colored nodes are active with u < 
0. The size of the node is proportional to its eigenvalue centrality. 
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FIG. 3: Robustness of quantized patterns. (A) Energy of the resulting quantized state with respect to 
the input signal amplitude u at the nearest and next-nearest neighbors of the best connected node in the 
system. When the amplitude is very small, the initial perturbation relaxes back to the trivial solution and 
no quantized state is formed (/). As the amplitude of the input signal is increased, fragile quantized states 
are formed (//). When the amplitude of the input signal is larger than a threshold value, a very robust 
quantized state is formed {Hi). Further increases in the input signal amplitude lead to the same quantized 
state, (insets) Visualization of the input signal in our network topology (the amplitude increases from 
left to right) as well as the resulting equilibrium state. (B) The energy of the resulting quantized state 
with respect to the ratio between the signal amplitude and the noise amplitude. Starting from the step¬ 
like input signal that gives the robust quantized state, we add random noise at the already perturbed 
neighborhood and we compute the energy of the resulting quantized state over 100 realizations. We use 
a Barabasi-Albert scale-free network of size V=200 and mean degree {k)=A. 
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FIG. 4: Global Turing patterns. Global patterns are possible when the non-active stationary solution is 
perturbed outside the parameter region of localized patterns (p < 0). The initial exponential growth of 
the perturbation is followed by a nonlinear process leading to the formation of stationary Turing patterns. 
(A) (left) The activation profile as a function of the node index i of a global stationary Turing pattern 
from direct simulation (blue crosses) is compared with the mean-field bifurcafion diagram. Black curves 
indicafe sfable branches while grey curves correspond fo unsfable branches of a single acfivafor-inhibifor 
sysfem coupled fo fhe compufed global mean fields. We sorf fhe node index i in increasing connecfivify k. 
Nodes wifh fhe same degree are sorfed wifh increasing fwo-jump connecfivify (see Insef). We 
use fhe same Barabasi-Albert nefwork model as in Fig. 2 and we sef fhe bifurcafion paramefer equal 
fo /X = —1/4. We have confirmed fhaf similar resulfs hold for larger nefwork sizes, (righf) Visualizafion 
of fhe global acfivify paffern on fhe nefwork topology. (B) The acfivafion profile as a funclion of fhe node 
index i of global sfafionary Turing pafferns from direcf simulation for bifurcation paramefer p = —0.25 
on an Erdos-Renyi random nefwork wifh size N = 1000 and mean degree {k) = 4 (blue curve) along 
wifh fhe sfable branch of fhe mean field approximation (black curve). (C) The acfivafion profile as a 
function of fhe node index i of global sfafionary Turing pafferns from direcf simulation on a Barabasi- 
Alberf scale free nefwork wifh fhe same mean degree and fhe same number of node as fhe Erdos-Renyi 
of B. We sorf fhe node index in increasing connecfivify k and fwo-jump connecfivify k^'^\ 
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